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PREFACE 


The  purpose  of  this  work  was  to  obtain  the  solution  to  three 
supersonic  flow  problems  using  three  different  numerical 
techniques . 

First,  a  shock/boundary  layer  problem  is  solved  using 
MacCormack ' s  explicit  technique.  Then,  using  the  same  technique 
a  shrouded  rocket  nozzle  problem  is  solved.  These  two  problems 
showed  that  the  explicit  scheme  required  many  minutes  of  computer 
time  to  solve. 

In  order  to  explore  more  efficient  codes  to  solve  the  shroud 
problem,  space  marching  algorithms  were  studied.  A  space 
marching  algorithm  using  flux-splitting  in  the  streamwise 
direction  was  applied  to  an  approximate  form  of  the  Navier-Stokes 
equation.  Flux-splitting  combined  with  a  global  Iteration 
approach  should  allow  the  shroud  problem  to  be  solved  with  a 
space  marching  algorithm.  The  flux-splitting  code  was  applied  to 
t  -to  supersonic  flow  problems  and  very  good  results  were  obtained. 

This  work,  of  course,  is  not  the  result  of  m£  labor  alone. 

Dr  Elrod  helped  me  a  great  deal  on  the  shroud  problem.  Major 
Hodge  showed  me  how  to  use  the  explcit  code  and  familiarized  me 
with  the  Cray  XMP.  Dr  Shang,  of  the  Flight  Dynamics  Lab,  has 
helped  me  understand  the  nuiances  of  MacCormack's  explicit 
scheme.  Finally,  I  have  to  thank  my  wife,  Jackie,  for  having  the 
patience  to  endure  my  absence  from  her  for  these  past  few  months. 
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Abstract 


In  the  present  study,  solutions  to  supersonic  fluid  dynamic 
problems  using  different  numerical  techniques  are  given. 

First,  the  full  Navier-Stokes  equations  are  implemented  in 
an  explicit  scheme  to  solve  a  shock/boundary  layer  interaction 
problem  and  a  shrouded  rocket  nozzle  problem.  In  the  shroud 
problem,  Mach  contours,  pressure  contours,  velocity  vectors, 
shroud  wall  pressure  and  base  wall  pressure  were  plotted. 
Comparison  of  the  base  wall  pressure  with  experimental  values 
were  inconclusive.  Shroud  thrust  increment  calculations  showed 
that  the  shroud  only  produces  thrust  gain  at  ambient  pressures 
very  close  to  zero.  It  was  also  found  that  a  great  deal  of 
computer  time  was  needed  to  solve  these  problems  with  the 
explicit  method. 

To  arrive  at  a  more  efficient  scheme,  an  approximation  was 
made  to  the  full  equations.  The  approximation  is  valid  for 
supersonic  flows  with  thin  shear  layers  and  also  enables  a  more 
efficient  scheme  to  be  applied  to  the  equations. 

These  schemes  are  called  space  marching  schemes  and  allow  a 
solution  of  a  flow  without  separation  in  a  single  sweep  of  the 
flowfield.  Problems  with  separation  regions  can  also  be  solved, 
but  a  global  iteration  process  which  requires  multiple  sweeps  is 
needed. 

Two  space  marching  techniques  were  explored.  The  first, 
credited  to  Vlgneron,  splits  the  pressure  from  the  streamwise 
flux  vector.  This  code  was  written  for  single  sweep  solution  of 
a  flowfield  and  it  was  used  to  solve  a  subsonic  Couette  flow  and 


a  supersonic  boundary  layer  problem.  The  second  code,  which  is 
original  to  this  work,  uses  flux-splitting  on  the  streamwise  flux 
vector  and  implements  the  global  iteration  technique.  The  same 
two  problems  solved  using  Vigneror. *s  technique  are  solved  using 
the  flux  splitting  technique  and  the  shock/boundary  layer 
interaction  problem  is  also  solved. 

It  was  found  that  Vlgneron*s  technique  is  more  efficient  for 
solving  flows  without  separation.  The  flux-splitting  code, 
however,  was  found  to  be  about  fourteen  times  faster  than  the 
explicit  code  when  applied  to  the  shock/boundary  layer  interact¬ 
ion  problem. 

An  efficiency  parameter  in  the  form  of  CPU  time  divided  by 
the  number  of  grid  points  divided  by  the  number  of  iterations  was 
calculated  for  each  code.  For  Vigneron's  method  the  value  was 
4.7  x  10~4  r  for  the  explicit  method  the  value  was  2.2  x  10"*  and 
for  the  current  method  the  value  was  4.6  x  10-4 


I.  INTRODUCTION 


Engineers  are  continually  looking  for  ways  to  increase  the 
efficiency  of  rocket  engines  since  for  every  additional  pound  of 
thrust  gained  another  pound  of  payload  can  be  added  or  a  higher 
orbit  can  be  achieved.  One  way  of  increasing  the  thrust  is  to 
increase  the  effectiveness  of  the  nozzle.  Nozzles  are  designed  for 
a  specific  altitude  and  at  any  other  altitude  the  nozzle  thrust 
is  not  optimum.  Deployable  nozzles  are  one  method  of  increasing 
the  thrust  above  the  design  altitude  and  shrouding  the  nozzle  is 
possibly  another  method. 

Figures  1  and  2  show  the  basic  features  of  shrouded  rocket 
nozzles.  When  the  flow  from  the  nozzle (s)  is  overexpanded  the 
shroud  is  ineffective  since  the  exhaust  plume  will  not  attach  to 
the  shroud.  But,  when  the  flow  from  the  nozzle (s)  is 
underexpanded  the  exhaust  plume  will  expand  out  and  attach  to  the 
shroud  wall  essentially  providing  a  larger  exhaust  exit  area. 

There  is  also  the  possibility  that  the  shroud  will  be  short 
enough  so  that  the  expanded  flow  will  not  attach  to  the  shroud, 
but  that  case  will  not  be  studied  here.  Allowing  the  exhaust 
plume  to  further  expand  so  that  the  exit  pressure  is  closer  to 
ambient  pressure  means  that  the  exhaust  flow  is  nearly  optimum 
and  hopefully  the  thrust  is  a  maximum.  However,  as  will  be  seen 
in  Chapter  II,  many  factors  affect  the  thrust  of  the  shroud  and 


Shroud 


uinnninnnvi  in  uvnnwuvwwuiwiwu  uuuuw  uuuiuu  wuu1 

P 

$ 


Figure  1  Shrouded  Nozzle:  Overexpanded 
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the  addition  of  the  shroud  may  cause  a  net  thrust  gain  or  a  net 
thrust  loss  compared  to  the  unshrouded  nozzle. 

Goethert  (7)  was  the  first  to  study  the  effects  of  shrouded 
rocket  nozzles.  Actually,  Goethert  used  the  shroud  to  eliminate 
the  problem  of  nozzle  burn  through  which  occured  on  the  clustered 
nozzles  of  the  early  Polaris  missies.  After  Goethert,  Holmes  and 
Matz  (8)  studied  the  effects  of  shroud  shape,  shroud  size  and 
nozzle  spacing  on  the  shroud  flow  field.  More  recently,  Moran  (11) 
has  completed  experimental  studies  with  various  shroud  shapes  and 
nozzle  configurations  on  both  two  and  three  dimensional  nozzles. 

On  the  numerical  side,  Roache  and  Mueller  (13)  studied 
incompressible  and  compressible  flow  over  backsteps,  which  is 
similar  to  a  single  shrouded  nozzle.  Recently,  Bardina  and 
Lombard  (3)  have  applied  a  three-dimensional  split  coefficient 
implicit  code  to  shrouded  mult i -nozzles  and  successfully  predict¬ 
ed  the  characteristic  features  of  the  flow. 

The  primary  objective  of  the  present  work  is  to  develop 
numerical  algorithms  capable  of  efficiently  simulating  predomin¬ 
ately  supersonic  flowfields,  such  as  shrouded  rocket  nozzles. 
Numerical  simulation  of  the  flow  inside  the  shroud  will  normally 
give  more  detailed  information  than  data  obtained  from  experi¬ 
ment.  This  information  may  lead  to  an  accurate  understanding  of 
the  flow  phenomena  and  consequently  the  optimum  performance  of 
the  shroud  can  be  predicted. 

The  flow  problem  inside  the  shroud  of  the  rocket  nozzle  is 
formulated  using  the  complete  unsteady  Navier-Stokes  (NS) 
equations  in  terms  of  the  conserved  variables.  The  numerical 
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scheme  used  to  solve  this  problem  is  the  explicit  MacCormack 
method  (2:479-489). 

Results  are  first  obtained  for  one  special  case  for  which 
published  results  are  available;  this  consists  of  the 
shock/boundary  layer  interaction  problem  (4,18).  Comparison  of 
the  present  results  with  the  available  data  provides  a  test  of 
the  validity  and  correctness  of  the  present  analysis  and  solution 
procedure . 

Solutions  are  then  obtained  for  the  case  of  primary  interest 
in  the  present  study,  namely  the  flow  for  a  single  shrouded 
nozzle.  This  model  problem  corresponds  to  experiment  3B  from 
reference  11.  The  geometric  details  of  the  problem  are  given  in 
Appendix  A. 

Upon  completion  of  the  shroud  problem  it  was  found  that  the 
CPU  time  required  to  obtain  the  solution  was  excessive.  This 
results  from  the  nature  of  the  explicit  scheme  which  has  a  limit¬ 
ation  on  the  time  step  in  order  to  insure  a  stable  scheme. 
Therefore,  many  iterations  and  large  amounts  of  computer  time 
were  required  to  reach  a  steady-state  solution.  To  remove  the 
time  step  restriction,  fully  implicit  methods  have  been 
investigated.  The  implicit  methods,  however,  still  require  many 
iterations  to  reach  the  steady  state  and  consequently,  still 
require  large  computational  costs. 

In  an  effort  to  decrease  the  computational  costs  associated 
with  the  implicit  algorithms,  space  marching  procedures  have  been 
studied.  Partial  differential  equations  are  classified  according 
to  three  types  —  elliptic,  parabolic  and  hyperbolic.  In 
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elliptic  type  problems  the  solution  at  one  point  depends  on  the 
solution  at  all  other  points  in  the  domain.  The  problem  must 
then  be  solved  within  a  closed  domain  so  that  the  influence  of 
all  points  in  the  domain  is  taken  into  account.  In  parabolic  and 
hyperbolic  type  problems,  however,  the  solution  at  one  point  only 
depends  on  the  points  in  a  reqion  of  the  domain.  These  types  of 
problems  are  called  marching  problems  since  the  solution  advances 
outward  from  known  initial  conditions. 


Figure  3  Solution  domain  for  a  marching  problem 
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In  the  present  efforts,  two  types  of  marching  algorithms  are 
investigated  to  solve  the  steady,  approximate  Navier-Stokes 
equations  (AMS).  The  ANS  equations  differ  from  the  NS  equations 
in  neglecting  the  streamwise  diffusion  terms  and  the  time  depen¬ 
dent  terms.  Now,  when  the  flow  is  supersonic  the  ANS  equations 
are  parabollc/hyberbolic  in  nature,  which  allows  the  solution  to 
be  marched  in  the  streamwise  direction.  Unlike  the  parabolic 
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boundary  layer  equations,  the  AMS  equations  retain  the  interact¬ 
ion  between  the  viscous  and  inviscid  portions  o£  the  £low  field 
which  allows  the  AMS  equations  to  model  elliptic  flows. 

Since  the  ANS  equations  are  elliptic  when  the  flow  is 
subsonic,  then  some  method  must  be  found  to  suppress  the 
ellipticity  so  that  the  equations  may  be  marched.  Rubin  and  Lin 
(15),  Schiff  and  Steger  (16)  and  Vigneron,  Rakich  and  Tannehill 
(19)  have  successfully  applied  such  methods.  In  all  of  these 
methods,  the  streamwise  flux  vector  was  split  in  a  way  that 
treats  the  pressure  term  separately.  For  example,  in  the  Schiff 
and  Steger  method,  the  pressure  was  assumed  to  be  constant  in  the 
viscous  subsonic  layer  and  equal  to  the  value  from  the  adjacent 
supersonic  flow,  while  Rubin  and  Lin  indicate  that  setting  the 
streamwise  derivative  of  the  pressure  to  zero  can  be  a  valid 
approximation  for  flows  with  cold  walls.  In  any  case,  the  error 
is  confined  to  only  a  thin  sublayer  adjacent  to  the  wall. 
Vigneron,  Rakich  and  Tannehill  approximated  the  streamwise 
derivative  term  with  a  weighting  between  implicit  and  explicit 
differencing  that  depends  on  the  local  Mach  number. 

The  method  of  Vigneron,  Rakich  and  Tannehill,  known  as 
Vigneron' s  method,  was  used  to  solve  two  problems.  The  first 
problem  is  the  Couette  flow  inside  a  channel  and  the  second 
problem  is  the  supersonic  laminar  flow  over  a  flat  plate.  Be¬ 
cause  the  Couette  problem  is  a  fully  developed  flow,  the  stream- 
wise  derivatives  are  zero  and  the  ANS  equations  are  the  exact 
equations  for  the  problem.  Also,  the  Couette  flow  has  an  exact 
analytical  solution.  The  supersonic  boundary  layer  flow  has  been 


solved  by  Lawerence  and  Tannehill  (10)  and  provides  a  more  severe 
check  £or  the  code. 

A  new  approach  to  spatially  marching  the  ANS  equations  is 
implemented  in  the  present  work.  In  this  method,  referred  to  as 
the  current  method,  the  Van  Leer  flux-splitting  technique  (1)  is 
used  to  decompose  the  streamwise  flux  vector  into  two  flux 
vectors  which  separately  model  the  information  traveling  upstream 
and  the  information  traveling  downstream.  Such  splitting  would 
appear  to  be  better  than  the  pressure  splitting  method  since  it 
is  more  representative  of  the  physics  of  the  flow  field  and  would 
assure  that  the  upstream  traveling  information  does  not 
destabilize  the  marching  algorithm.  Solutions  are  obtained  for 
the  same  two  problems  solved  using  Vlgneron's  technique  and 
both  are  compared.  Also,  the  solution  to  the  shock/boundary 
layer  interaction  problem  is  obtained  and  compared  to  the  solut¬ 
ion  from  the  explicit  code  and  the  results  obtained  by 
Thomas  ( 18 ) . 

Finally,  a  qualitative  comparison  between  all  the  methods 
is  given. 
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II.  THEORY 


In  this  chapter  the  equations  governing  the  flow  problems 
are  examined.  For  the  shroud  problem  a  simple  control  volume 
analysis  will  reveal  the  important  properties  that  affect  the 
thrust  gain  of  the  shroud/nozzle  configuration.  Next,  the 
equations  to  be  used  in  the  numerical  schemes  are  presented; 
starting  with  the  general  form,  advancing  to  the  vector  form  of 
the  unsteady  Navier-Stokes  equations  and  finally  arriving  at  the 
approximate  Navier-Stokes  equations.  Then,  the  flux-splitting 
technique  used  in  the  current  method  is  discussed.  Finally,  the 
chapter  will  conclude  by  presenting  the  analytical  solution  to 
the  Couette  flow  problem. 


SHROUD  ANALYSIS 

Many  Important  features  of  this  problem  can  be  analyzed  by  a 
simple  control  volume  analysis.  Equations  for  the  thrust  of  the 
shroud  and  the  nozzle  are  derived  by  using  the  x-momentum  equation 


F,  =  I  u*  dm 


Assuming  that  the  flow  properties  and  ambient  conditions  are  a 
constant,  then  Eq. 1  can  be  applied  to  the  control  volume  (CV) 
shown  in  figure  4  to  give 

T,  =  mu,  -  (p,  -  p.)  A, 
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This  is  the  familiar  rocket  thrust  equation.  Now,  by  applying 
Eq.  1  to  the  CV  shown  in  figure  5,  another  equation  for  the  thrust 
of  the  shroud  can  be  obtained 


T.  =  muB  -  (pb  -  pjAb  -  (pB  -  pjAn  -  f  (3) 

Recognizing  that  the  thrust  of  the  nozzle  alone  is  given  by 

T„  =  muB  *  (pB  -  p.)  Ab  ( 4 ) 


then  Eq  3  can  be  rewritten  as 


AT  =  T.-T.  =  (pk  -  pj  Ab  -  f 


(5) 


Eq  5  shows  the  Importance  of  the  shroud  base  pressure  in  in¬ 
creasing  the  thrust  of  the  nozzle.  If  the  shroud  base  pressure 
force  is  higher  than  the  ambient  pressure  force  plus  the  shear 
force,  then  the  thrust  increment  is  positive;  otherwise  the 
thrust  Increment  is  zero  or  negative. 

GOVERNING  EQUATIONS 

The  integral  analysis  given  above  can  only  give  Information 
on  the  gross  properties  of  the  flow,  such  as  the  thrust.  For  a 
more  detailed  analysis,  however,  the  equations  of  motion 
expressing  the  conservation  of  mass,  momentum  and  energy  are 
required  in  the  differential  form.  In  generalized  coordinates 


the  conservation  o£  mass  equation  Is 


dp/d  t  V.  (/>V)  =  o 

The  conservation  o£  linear  momentum  equation  is 


(6) 


I 

i 

t 

t 


i 
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p(DV/Dt)  =  -Vp  +  V*  r  +  pB 


The  conservation  o£  energy  equation  is 

p(De/Dt)  =  -p(V  (8) 


where  ♦  is  the  dissipation  function  and  is  given  by 

*  =  T  :  (V  V) 


In  order  to  complete  this  set  of  equations  the  pressure, 
the  stress  tensor  and  the  heat  conduction  vector  must  be  related 
to  the  other  variables.  Only  perfect  gases  will  be  considered 
here  and  thus  the  following  equation  of  state  will  be  used 


p  =  pRT 


(10) 


where  R  is  the  universal  gas  constant.  Since  only  Newtonian 
fluids  will  be  studied,  Stoke's  law  of  friction  can  be  used 
to  relate  the  shear  stress  tensor  to  the  gradient  of  the  velocity 
vector  as 
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r  =  X V  •  V  I  /i[W  +  (VV)T] 


(11) 


Finally/  the  heat  flux  vector  is  related  to  the  gradient  of  the 
temperature  through  Fourier's  law  of  heat  conduction 

q  =  -  kVT  (12) 

where  the  conductivity/  k,  is  a  function  of  the  temperature. 

For  numerical  solution  these  equations  can  be  expressed  in 
vector  form  in  cartesian  coordinates.  Also,  to  simplify  the 
finite  difference  algorithms  for  problems  requiring  complex 
grids,  the  equations  are  written  in  terms  of  a  computational 
space.  Eq  6  through  Eq  12  are  then  written  as 

a/at(u/j)  -  -  fcFrPJD 


(13) 


T„  =  (2/3)ix[2({.U^7J,UT?)-({yV^tJ,Y7J)] 

T„  =  (2/3M2UIv^^,vf|)-{^ur,,*un)l 


qf  =  -k  (tsT|^i|sT1|) 
q,  =  -k 


For  most  high  Reynold's  number  flows  solved  by  finite 
differences,  an  approximate  form  of  the  Navler-Stokes  equations 
(AMS)  can  be  solved.  The  ANS  equations  to  be  used  are  sometimes 
incorrectly  called  the  "parabolized"  Navler-Stokes  equations 
although  they  are  actually  parabolic/hyperbolic.  The  ANS  equat¬ 
ions  are  derived  from  the  NS  equations  with  the  following  reason¬ 
ing.  The  grid  densities  allowed  by  computer  memory  and  cost 
limitations  usually  only  allow  the  viscous  terms  to  be  accurately 
resolved  in  one  coordinate  direction.  The  viscous  terms  in  the 
direction  normal  to  a  surface  are  needed  for  such  quantities  as 
shear  forces  and  energy  transfer  rates.  The  streamwlse  viscous 
terms,  then,  are  not  accurately  resolved  even  when  the  NS  equat¬ 
ions  are  solved.  The  ANS  equations  are  arrived  at  by  simply 
neglecting  the  streamwlse  viscous  terms  and  the  time  dependent 
terms  of  the  full  Navler-Stokes  equations  written  in  Cartesian 
coordinates  and  can  be  expressed  as 


dEi’/dx*  -  dFi’/d/  -  dFf*/dy*  =  0 


E,*  = 


t  I 

P  u 

*  *3 
/>  U  ■* 

*  *  a 
P  U  Y 


Fi 


fi  U  V 
/>*¥**  D* 

(Ft*  -  pV 


F?*  =  /i*/ReL 


0 

du*/dy* 

4/3  dv*/dy* 

u*du*/dy*  (4/3)v*dv*/dy* 


[(7-l)M,*Prr'  aT*/dy‘ 


Where  the  Reynold's  number  Is  defined  by 


Ret  =  P|U,L/p, 


A/ 


l 


and  the  equation  of  state  is  now 

9 

t  =  />*t7km,j 


i  A 


,v 

.V 


H 


The  variables  are  non-dlmensionalized  as  follows 

u*  =  u/u,  v*  =  v/uf  x*  =  x/L 

y*  =  y/L  p*  =  p/p  i  e  =  e/u»* 

(28) 

P*  =  p/p|Uf2  $i  =  WMt  t*  =  T/T» 
a*  =  a/uj 

In  many  problems  the  viscous/inviscid  regions  Interact  and 
to  capture  this  interaction  with  the  fewest  number  of  grid 
points,  stretching  must  be  applied  in  the  normal  direction.  To 
facilitate  this  stretching,  the  equations  are  written  in 
computational  coordinates  with  only  the  normal  coordinate 
transformed.  The  ANS  equations  can  be  easily  modified  for  these 
coordinates  as  follows.  The  normal  direction  in  computational 
coordinates  will  be  designated  as  eta.  The  normal  direction 
derivative  in  physical  space  can  then  be  written  in  terms  of  the 
normal  direction  derivative  in  computational  space 

djdy*  =  dri/dy*  d/drj 
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substit luting  into  Eq  22  yields 


dEiVdx*  -  =  0 

dividing  by  f}yV  gives 


(1/I|f.)  dR*/dx*  +  d?'/d7)  -  dF,Vdt|  =  0 


where 


PT*  =  ti\,/ Re^ 


0 

du*/di ? 

(4/3)dyVat| 

udu/frn  ♦  (4/3)v'd//di>  •*-  [(7-l)M,*Prr‘  dT*/dij 


(30) 


These  equations  are  used  in  the  current  method. 


Flux  Splitting 

Since  the  streamwise  flux  vector  of  the  ANS  equations  only 
contains  invlscid  terms,  then  this  vector  may  be  "split"  into  two 
flux  vectors.  One  flux  vector,  called  the  plus  flux  vector, 
models  information  propagation  downstream.  The  negative  flux 
vector  models  information  propagation  upstream. 

Starting  with  the  derivative  of  the  streamwise  flux  vector 

dtfjdx* 
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This  can  be  written  as 


[q]  au*/a*' 


where 


[Q]  «  dBiVau* 


Now  gj*  has  the  property  (2:282) 


e;  =  rau 


* 


(31) 


(32) 


A  similarity  transformation  also  exits  so  that  (2:281) 

M  [Q1  [R]  =  M  (33) 

where 

[L]  =  matrix  of  left  eigenvectors, 

[RJ  =  matrix  of  right  eigenvectors, 

[A]  =  matrix  of  eigenvalues 

The  eigenvalue  matrix  can  be  split  into  two  parts.  One  part 
contains  the  positive  eigenvalues  and  the  other  part  contains  the 
negative  eigenvalues.  This  is  written  as 

M  =  [*r~w  ,34, 

Substituting  Eq  34  into  Eq  33 
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w  [qj  [Rj  =  [xr  -  [xr 


Solving  for  IQ] 


mi  =  iLr‘  [xr  [Rr1  -  dr1  [xr  [Rr 
=  for  -  tor 


The  streamwise  flux  vector  can  then  be  written 


Ei*  -  wr  u*  -  wr  u* 

=  Ei*"  -  Ej*~ 


(36) 


The  exact  elements  of  E^*  and  Ej~*  may  differ.  The  form  used  in 
this  work  was  derived  by  Van  Leer  and  is  given  by 


E^«.v* 

E^liM.a(7-l)u*-2»,r/[2(7,-l)]  -  y**/2} 


where 


p\1(l/2)  (M.-lff 


(38) 


vxtx 


When  the  streamwlse  Mach  number  ie  greater  than  or  equal  to 
one  then 

B,**  =  B,* 

and  when  the  streamwlse  Mach  number  is  less  than  or  equal  to 
negative  one  then 
-*• 

When  combined  with  the  proper  finite  difference 
approximations,  flux  splitting  provides  a  realistic  model  of  the 
physics  of  a  flow.  Ej**  models  the  information  propagation  in  the 
downstream  direction  and  Ej”*  models  the  information  propagation  in 
the  upstream  direction.  Appropriate  finite  differences,  then 
will  take  the  correct  Information  propagation  into  account  when 
approximating  the  derivatives  ofBj**and  Ej~*  •  For  the-Ej** 
derivative  a  backward  difference  would  be  used  and  for  the  Ej"* 
derivative  a  forward  difference  would  be  used.  An  example  of  a 
backward  difference  is 

(du/dx)jj  -  [Ujj-Uj_,J/AX 

and  an  example  of  a  forward  difference  is 

(du/dx)y  ~  [Uj^y-UjJ/AX 

This  splitting  method  will  be  applied  to  the  current  method 
described  in  Chapter  III. 


Couette  Flow  solution 

A  Couette  flow  la  an  Incompressible  flow  which  consists  of  a 
stationary  lower  plate  and  a  moving  upper  plate  with  a  fluid 
between  the  plates.  For  the  purposes  of  this  work  the  fluid  is 
air  and  the  flow  is  considered  two-dimensional  and  fully 
developed.  With  these  assumptions  the  x-momentum  equation 
of  the  NS  or  ANS  equations  reduces  to 

dV/dy**  =  0 

With  the  boundary  conditions 

u  V=0)  =  0  ;  u*(y*=l)  *  1 

the  solution  is 


u  V)  = 


t 


y 


(39) 


The  energy  equation  is 


(du'/dy*)1  -  (d*T*/dy**)  =  0 


and  the  boundary  condtions  are 

T*(y*=0)  =  1  ;  T*(y*=l)  =  1 
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The  important  concepts  of  solving  the  governing  equations 
using  finite  difference  methods  are  discussed  in  this  chapter. 
First#  the  governing  equations  are  examined  carefully  with  the 
aim  of  selecting  an  efficient  numerical  scheme  for  their 
solution.  Thereafter,  the  finite  difference  form  of  the 
equations  is  presented,  together  with  the  initialization  used  for 
all  variables.  The  procedure  for  solving  the  resulting  algebraic 
simultaneous  equations  is  described  and  some  details  about  the 
implementation  of  the  boundary  conditions  are  discussed. 

Inherently,  the  governing  equations  are  an  extremely  stiff, 
non-linear  system  of  equations  and  therefore  their  numerical 
solution  deserves  special  care.  The  most  widespread  method  of 
solution  is  the  use  of  finite  difference  discretization.  Simply 
put,  the  finite  difference  method  converts  partial  differential 
equations  (PDE)  into  a  system  of  algebraic  equations.  These 
equations  are  then  written  at  the  discretized  grid  points  used  to 
approximate  the  flow  domain.  The  system  of  equations  can  then  be 
solved  explicitly  or  implicitly. 

In  an  explicit  scheme  the  finite  difference  form  of  the 
differential  equation  is  written  so  that  only  one  unknown  appears 
in  the  equation.  For  example,  the  heat  equation 


du/dt  =  a  dVdx* 


(41) 
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An  explicit  formulation  is 


I 


(Uj*“  -  Uj*)/At  =  (a/(Aif)  [„w*  -  2U(*  *  uH*] 


9 

> 


% 


ss 


where  n  is  the  time  index  and  j  is  the  space  index.  Because  an 
initial  condition  must  be  specified,  then  the  values  at  n  will 
always  be  known.  Only  one  unknown,  Uj*-*-1  ,  appears  in  the  equation. 
Explicit  methods  are  generally  easy  to  code,  but  require  very 
small  time  steps  to  satisfy  stability  restrictions.  Therefore, 
many  iterations  and  large  computer  times  are  required  to  reach  a 
steady-state  solution.  This  can  be  a  severe  problem 
when  very  fine  grids  must  be  used  to  capture  flow  details. 

An  Implicit  finite  difference  equation  for  the  heat  equation 
could  be  written  by  evaluating  the  terms  on  the  right  hand  side 
of  Eq  42  at  the  n+1  time  level.  The  implicit  schemes,  then 
require  the  simultaneous  solution  of  several  equations.  For  a 
linear  system  of  equations  implicit  methods  are  unconditionally 
stable.  However,  to  be  able  to  solve  the  discretized  form  of  the 
governing  equations  implicitly,  they  must  first  be  linearized. 
This  leads  to  a  matrix  system  of  equations  which  can  be  solved  to 
obtain  the  variables  at  each  grid  point,  but  the  method  is  no 
longer  unconditionally  stable.  The  maximum  time  step  is  much 
larger  than  allowed  for  the  explicit  methods  so  that  fewer 
iterations  and  lower  computer  times  are  required.  Also,  the 
explicit  methods  may  fail  for  problems  having  very  steep  grad¬ 
ients  so  the  only  choice  is  to  use  an  implicit  method. 


« * 
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li  SSBllglt  Method 


The  explicit  method  used  in  this  work  is  a  two-dimensional 
version  of  MacCormack's  predictor-corrector  scheme ( 2 : 479-489 ) . 
This  scheme  was  chosen  because  it  was  readily  available  and  had 
been  successfully  used  in  earlier  AFIT  studies(14).  The 
empirical  stability  constraint  on  this  scheme  is  given  in 
cartesian  coordinates  by  (2:484) 


At<  -  2/Rsa  ) 


(43) 


where  sigma  is  approximately  0.9  and  (AtJon,  is  the  lnvlscld  CFL 
condition  given  by  (2:484) 


(AtJow.  <_  [M/x  *  M/y  ♦  <V(ax)*  +  i/(Aj)7/*r‘ 


(44) 


The  particular  code  used  in  this  study  was  obtained 

from  the  Computational  Aerodynamics  Group,  Flight  Dynamics  Lab, 

Wr ight-Patterson  Air  Force  Base.  The  code  has  been  vectorized 
for  the  Cray  XMP  supercomputer  and  contains  a  Baldwln-Lomax 
arithmetic  turbulence  model,  although  all  the  problems  considered 
here  will  be  laminar.  The  explicit  code  discretizes  the  dimen¬ 
sional  form  of  the  uns4"  ady  NS  equations  in  computational 
coordinates  with  all  variables  in  English  units.  For  more 
details  on  MacCormack’s  explicit  scheme  see  reference  2  and  for 
the  details  on  the  code  see  reference  14. 


.V  •  j  v  V  ’  V  V  /  -*.v.  /. 


•  »< 


% 


M 

i 


I 


rj 


V 


/  ' 


a 


.V 


Space  Marching  Algorithms 


The  first  method  follows  the  method  used  by  Vlgneron  to 
solve  the  AMS  equations  using  a  space  marching  scheme.  This 
allows  the  solution  of  a  predominately  supersonic  flow  to  be 
obtained  with  a  single  sweep  of  the  problem  domain.  In  the 
subsonic  portions  of  the  domain  information  traveling  upstream 
must  be  suppressed,  and  Vlgneron  achieved  this  by  "splitting"  the 
streamwise  pressure  gradient  from  the  streamwise  flux  vector.  Eq 
22  would  then  be 


dE^/dx*  -  aP'/dx*  -  d?i!dy  -  dFt7dy'  =  0 


(45) 


where 


5*-  = 


•  • 

P  u 

•  ••  • 
P  u  ctp 
•  •  • 

PUT 


(K,*  *  py 


(46) 


and 


P*  = 


0 

0 


(47) 


Vlgneron's  stablillty  analysis  showed  that  for  the  solution 
to  be  marched  in  the  streamwise  direction  omega  must  satisfy 


•  <  7M,7[1  -  (7-1JM,*] 


(48) 
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For  computational  purposes  a  safety  factor  is  applied  to  Eq  48 


4*  =  <T7M,*/ri  -  (7-  1JM/J 


(49) 


where  sigma  is  usually  0.9.  The  Beam-Warming  finite  difference 
equation  Is  given  by 


a'B  -  txdfdx  (a'B)  ♦  txd/dx  (B1) 


(50) 


where 


a'B  «  E1*1  -  B* 


(51) 


This  is  the  Euler  Implicit  form  which  is  first  order  accurate  in 
the  inarching  (streamwise)  direction  and  second  order  accurate  in 
the  normal  direction.  When  Eq  50  is  applied  to  Eq  45  the  stream- 
wise  pressure  gradient  terra  is  treated  as  an  explicit  terra  and 
taken  to  the  right  hand  aide  (RHS).  This  term  is  then  backward 
differenced  to  retain  the  space  marching  technique.  However, 
when  this  is  done  the  streamwise  step  size  must  meet  the 
stability  requirement 


where  p  is  the  wave  number.  This  requirement  can  be  a  severe 
restriction  for  low-speed  flows.  This  restriction  can  be  removed 
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using  two  methods  —  neglecting  the  streamwise  pressure  gradient 
term  and  specifying  the  pressure  field.  Neglecting  the  stream- 
wise  pressure  gradient  is  an  obvious  way  to  remove  the  restric¬ 
tion/  but  the  method  must  then  be  used  only  for  flows  with  small 
or  no  streamwise  pressure  gradient.  Specifying  the  pressure 
field  allows  the  streamwise  pressure  gradient  term  to  be  forward 
differenced  and  again  removes  the  stability  restriction. 

However/  if  the  pressure  field  is  not  exactly  known,  then  the 
method  must  be  modified  for  global  iterations  where  the  field  is 
swept  until  the  correct  solution  is  reached.  This  diminishes  the 
attractiveness  of  the  single-sweep  marching  technique,  but  only  a 
few  sweeps  are  usually  required.  The  code  developed  for  the 
present  work  using  Vigneron's  method  will  be  referred  to  as  PANS. 
Vigneron's  method  closely  follows  the  current  method  except  that 
the  PANS  code  is  a  single  sweep  code  while  the  current  method  is 
written  for  global  iterations. 

The  goal  here  is  to  develop  a  method  which  can  successfully 
be  applied  to  large  separation  regions.  Because  the  flux-split- 
ting  correctly  models  the  characteristic  propagation  of  informat¬ 
ion  in  both  subsonic  and  supersonic  flows  it  would  be  a  better 
candidate  for  a  successful  method.  However,  to  correctly  apply 
the  flux-splitting  a  forward  difference  must  be  applied  to  the 
negative  flux  term.  This  means  that  the  marching  scheme  must  be 
modified  for  global  iterations  since  the  forward  difference 
requires  data  that  has  not  been  calculated.  To  solve  this 
problem  the  negative  flux  term  will  be  calculated  using  "old" 
data,  that  is,  data  obtained  from  the  previous  solution  sweep. 
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The  derivation  of  the  current  algorithm  follows  that  given  for 
Vigneron's  method  in  reference  2.  Putting  Eq  22  into  delta  form  and 
substituting  into  Eq  50  gives 


A*Ei*  =  -Axld/dy*  (a'F/  -  a*Ft*)]  -  hx[d/dy*  (F,*1  -  F,n)] 


From  Eq  36  the  following  equation  can  be  written 

a'E/  *  A1  Bf*  -  A1^"* 


Substituting  Eq  36  into  Eq  53  results  in 


a'E,**  a'E,*-  =  -Ax*[d/dy*  (a'F,*  -  a‘F/)]  -  A*[d/dy*  [F*  -  F,' 


)] 


(54) 


This  equation  is  linearized  according  to 

a'f;-  [Rfi'U-  (55) 

a'F/~  [WJ's'U'  (56) 


where 


[R]1  *  (dF(Vau7 


(57) 


[Wf  *  (dFf7dU*y 


(58) 


According  to  the  flux-splitting  method  the  positive  streamwise  flux 
vector  for  subsonic  streamwise  Mach  number  is  written 

a'b“  ~  wrW  ,  s, , 


and  when  the  streamwise  mach  number  is  supersonic 

a'E’WE*  ~  [Q^a'U*  <60> 

The  elements  of  [Q] ,  [QF  ,  ( R 3  and  [V]  are  given  in  Appendix  B. 
Now  Eq  54  becomes 

[gr  w  -  a'e.*-  =  -«*[d/d/([R]W  -  [W?a'u*)]  (61) 

-A x\d/0,'  (F;  -  F/)] 


Taking  the  Implicit  terms  to  the  left-hand-side  (LHS)  and  the 
explicit  terms  to  the  right-hand-side  (RHS )  results  in 


[QTW  -  tx\d/dy*  ([R^U*  -  [W]VU*)] 

=  -Axld/d/  {?:  -  F,*)]  -  A V 


(62) 


The  first  partial  derivatives  are  approximated  by  second  order 
central  differences,  for  example 


°f°y aR]W) ~  . ([r]W)w}/(2a,-)  (63) 


and  for  the  second  derivative  terms 


mr  I M w  (V4'U0]  -  fa,-  * 

-  (M|’*m‘wX(wai01*J1  -  (va'U, %.,]/[!!(*/)») 


(64) 


This  gives  the  LHS  a  block  trl-dlagonal  structure  whose  elements  are 
given  In  Appendix  C.  The  negative  split-flux  vector  on  the  RHS 
is  evaluated  using  a  first  order  forward  difference 


bU%' 


(65) 


The  viscosity  coefficients  were  evaluated  using  Sutherland's  law 


n'  -  |*V  {T*  *'7(T*  *  T,*)} 


(66) 


The  current  code,  referred  to  as  FANS,  Is  carried  out  as  follows 
Step  1:  Specify  the  Initial  conditions  for  the  whole 
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flowfield.  For  the  PAMS  code  only  the 
conditions  at  1=1  need  to  be  specified. 

Step  2:  Calculate  Ej”*for  the  whole  flowfield.  Skip  this  step 
for  the  PANS  code. 

Step  3:  Calculate  the  streamwise  Mach  number  and  the  viscosity 
coefficient  for  each  J  (normal)  point  at  the  I 
station.  For  the  PANS  code  calculate 
omega  Instead  of  the  streamwise  Mach  number. 

Step  4:  Calculate  the  p*  vector  for  each  J  at  I 

Step  5:  Calculate  the  derivative  of  P*  with  respect  to  y  for 
each  J  at  I 

Step  6:  Calculate  the  derivative  of  p^*  with  respect  to  y  for 
each  J  at  I 

Step  7:  Calculate  the  RHS  vectors  for  each  J 

Step  8:  Calculate  the  linearization  Jacobians.  For  the  FANS 
code  if  the  local  streamwise  Mach  number  is 
less  than  one  then  calculate  [q r  instead  of  [Q] 

Step  9:  Calculate  the  LHS  matrices 

Step  10:  Solve  for  delta-U  at  1+1  using  the  block  tri- 
diagonal  matrix  solver 

Step  11:  Update  the  solution  using 

U*w  =  y»i  A‘u* 

Step  12:  Calculate  the  explicit  boundary  conditions 

Step  13:  Repeat  Steps  3  through  12  until  the  solution  for  the 


last  streamwise  station  is  obtained.  The  first  sweep 
of  the  flowfield  is  finished  and  the  PANS 


code  stops  here 

Step  14:  Go  to  Step  2  and  repeat  the  process  until  the 
solution  has  converged. 

Please  note  that  the  PAMS  code  is  not  written  in 
computational  coordinates.  Therefore,  a  grid  does  not  need  to  be 
generated;  simply  specify  the  streamwise  step  size,  the  number  of 
points  in  the  streamwise  direction,  the  normal  step  size  and  the 
number  of  points  in  the  normal  directibn.  The  FANS  code  does, 
however,  require  a  grid  to  be  generated  for  the  normal  coordin¬ 
ate  . 

Although  both  codes  were  run  on  the  Cray  XMP,  no  effort  was 
made  to  vectorize  the  codes.  In  fact,  vector ization  was  not 
implemented  during  the  compilation  of  the  codes. 


IV  Results  and  Discussion 

In  this  chapter  each  problem  that  was  solved  will  be  briefly 
discussed.  The  problems  are  the  following: 

-  Shock/Boundary  Layer  Interaction  (explicit) 

-  Single  Nozzle  in  a  Short  Shroud  (explicit). 

-  Couette  Flow  (PANS,  FANS) 

-  Supersonic  Boundary  Layer  (PANS,  FANS) 

-  Shock/Boundary  Layer  Interaction  (FANS) 

A  brief  description  of  each  problem  will  be  followed  by  the 
results  and  finally  a  discussion  of  the  signifigance  of  the 
results . 

Shock /Boundary  Laver  Interaction  (explicit ) 

This  problem  consists  of  a  shock  wave  impinging  on  the 
boundary  layer  over  a  flat  plate  as  shown  in  figure  6.  The 
results  obtained  here  will  be  compared  to  the  results  of  Thomas 
and  Walters(18)  and  also  the  results  of  Beam  and  Warming(4). 

The  computational  grid,  shown  in  figure  7,  is  the  same  as 
used  by  Beam  and  Warming  but  the  upper  boundary  was  Increased  to 
allow  the  leading  edge  shock  to  exit  through  the  downstream 
boundary  resulting  in  a  32x62  cartesian  grid.  The  boundary 
conditions  for  the  supersonic  inflow  were  freestream  up  to  the 
point  where  the  shock  entered  the  domain  and  then  the  inviscid 
shock  jump  conditions  from  that  point  to  the  top  of  the  domain. 
The  conditions  along  the  top  of  the  domain  were  also  set  to  the 
shock  jump  conditions  and  along  the  outflow  boundaries  the 
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conditions  were  extrapolated  from  values  of  Interior  points. 

Along  the  bottom  of  the  domain  symmetry  conditions  were  used 
before  the  leading  edge  of  the  plate  and  on  the  plate  no-slip 
and  constant  wall  temperature  conditions  were  applied.  The  wall 
temperature  was  set  to  the  free-stream  temperature  of  460 
Rankine . 

Piqure  8  shows  the  pressure  contours.  Figure  9  shows  the 
pressure  along  the  bottom  of  the  domain  and  Figure  10  shows  the 
friction  coefficient  along  the  plate.  The  results  agree  well 
with  the  solution  obtained  by  Beam  and  Warming  (4).  However,  a  few 
problems  should  be  pointed  out.  First,  a  large  amount  of  damping 
was  necessary  to  achieve  a  smooth  solution.  Next,  as  shown  by 
Thomas  and  Walters  (18),  a  grid  density  of  60x90  points  is  necessary 
to  accurately  resolve  the  separation  region.  The  grid  used  by 
the  explicit  code  wa3  32x62  points  and  required  22,000  iterations 
with  a  run  time  of  968  seconds  on  the  Cray  XMP.  The  efficiency 
parameter,  CPU  time  in  seconds  divided  by  the  number  of  grid 
points  divided  by  the  number  of  iterations,  was  2.2x10“*  .  Only 

one  streamwise  station  exhibited  the  negative  streamwise  velocity 
characteristic  of  the  separation  region.  Finally,  the  separation 
region  of  the  flow  was  very  slow  in  developing  which  reveals 
another  weakness  of  the  explicit  code  --  slow  convergence  of  the 
solution  in  subsonic  regions. 

Single  Nozzle  in  a  Short  Shroud 

The  problem  to  be  solved  now  is  the  two-dimensional  single 
nozzle  in  a  short  shroud,  experiment  3B  from  reference  11. 
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Figure  10  Shock/Boundary  Layer:  Friction  coefficient  using 

the  explicit  code 
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This  is  a  cold  air  flow  where  a  chamber  was  pressurized  to  13680 
pounds  per  square  foot  (psf)  and  the  exhaust  chamber  was 
evacuated  to  about  36  psf.  The  pressure  chamber  exhausted  into 
the  test  section  where  the  nozzle  block  and  shroud  block  were 
located.  The  exhaust  chamber  pressure  would  gradually  increase 
during  the  experiment  resulting  in  an  unsteady  flow  where  the 
nozzle  exhaust  plume  would  eventually  separate  from  the  shroud 
wall.  The  nozzle  had  a  4:1  exit  area  to  throat  area  with  an 
isentropic  exit  Mach  number  of  2.94.  More  detatils  on  the 
dimensions  are  given  in  Appendix  A.  For  this  numerical  solution 
the  exhaust  chamber  pressure  will  be  assumed  to  be  a  steady  36 
psf  so  that  a  steady  solution  with  the  exhaust  plume  attached  to 
the  shroud  wall  can  be  computed. 

In  order  to  keep  the  cost  and  computer  time  to  a  minimum  the 
number  of  grid  points  and  the  location  of  the  grid  points  in  the 
domain  must  be  carefully  selected.  The  flow  domain  was  chosen  to 
be  a  two-dimensional  cartesian  coordinate  system  with  x 
representing  the  streamwise  direction  and  y  the  normal  direction. 
The  reference  length  was  chosen  to  be  the  length  of  the  shroud  -- 
two  inches.  Because  of  the  large  separation  region  in  this  flow 
a  45x45  grid  was  chosen  as  the  best  trade-off  between  run  time 
and  accuracy.  To  improve  the  accuracy  of  the  results  the  grid 
points  were  clustered  around  high  gradient  regions  and  regions  of 
interest.  In  this  case,  the  high  gradient  region  was  near  the 
exit  of  the  nozzle  and  the  base  wall  was  the  region  of  interest 
since  this  is  where  the  pressures  from  the  experiment  were  taken. 
After  about  five  short  computer  runs  the  grid  in  Figure  11  was 


( 

t 

( 

I 

chosen  as  the  best  for  the  problem. 

Next,  the  boundary  conditions  will  be  discussed.  There  are 
five  important  regions  where  boundary  conditions  must  be  applied 
in  this  problem  --  the  shroud  wall,  the  base  wall,  the  nozzle 
exit,  the  symmetry  line  and  the  shroud  exit.  For  both  the  shroud 
wall  and  the  base  wall  no-slip  and  constant  wall  temperature  condit¬ 
ions  were  used;  to  calculate  the  density  at  the  wall  the 
normal  gradient  of  the  pressure  was  assumed  to  be  zero.  The  nozzle 
exit  (shroud  inlet)  conditions  were  taken  to  be  the  exit  condit¬ 
ions  for  an  isentropic  quasi-one-dimensional  nozzle  with  stagnat¬ 
ion  conditions  equal  to  the  conditions  in  the  pressure  chamber  of 
the  experiment.  In  order  to  Increase  the  grid  density  the  flow 
was  assumed  to  be  symmetrical  about  the  centerline  of  the  shroud 
so  that  the  upper  boundary  of  the  domain  is  a  symmetry  line. 

Symmety  line  conditions  were  set  to  be  zero  gradient  except  for 
the  normal  velocity  component  which  was  set  to  zero.  Finally, 
since  the  shroud  exit  is  a  predominately  supersonic  flow  then  the 
shroud  exit  conditions  can  be  calculated  from  interior  points  by 
simple  extrapolation. 

To  arrive  at  realistic  initial  conditions  the  problem 
must  be  thought  of  as  an  unsteady  problem  which  reaches  a  steady- 
state  solution.  The  initial  conditions  would  then  be  that  of 
the  exhaust  chamber  before  the  pressure  chamber  valve  was 
released,  which  were  taken  to  be  the  following: 

Po  =36  psf 
Uo  =  0.0  ft/sec 
▼o  = 


0.0  ft/sec 


460  Rankine 


l'^4  »j 
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T«  = 


Po  =  0.000045553  slugs/cu  ft 


except  at  the  nozzle  exit 


The  results  of  the  shroud  computations  are  presented  in 


Figures  12  through  16.  The  Mach  contour  plot  (Figure  12), 


pressure  contour  plot  (Figure  13)  and  velocity  vector  plot 


(Figure  14)  reveal  the  Important  features  of  the  flow.  The  Mach 


contour  plot  shows  the  expansion  of  the  nozzle  plume  from  a 


nozzle  exit  Mach  number  of  approximately  three  to  a  shroud  exit 


Mach  number  of  approximately  five  and  a  half.  Somewhat  harder  to 


see  is  the  compression  shock  formed  when  the  exhaust  plume 


strikes  the  shroud  wall.  This  feature  is  located  in  the  lower 


right  hand  corner  of  the  plot  and  shows  a  drop  in  Mach  number 


from  5.5  to  2.5  across  the  shock. 


Details  of  the  seperation  region  are  also  revealed  in  the 


Mach  contour  plot.  This  region  has  a  very  complex  flow  pattern 


containing  many  eddies  and  a  surprisingly  wide  range  of  Mach 


contours.  For  instance,  embedded  in  the  predominately  subsonic 


separation  region  is  a  supersonic  region.  This  region  stretches  from 


the  plume/wall  interaction  area  towards  the  shroud  base.  This 


phenomenon  also  shows  up  in  the  velocity  vector  plot.  Here,  the 


arrows  close  to  the  plume/shroud  interaction  region  are  seen  to 


point  towards  the  shroud  base.  Moving  closer  to  the  base,  the 


arrows  gradually  lengthen,  indicating  a  fluid  velocity  increase. 


until  the  flow  turns  parrallel  to  wall  and  the  velocity 


decreases . 


Both  the  velocity  vector  plot  and  the  pressure  contour  plot 
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where  w  Is  the  width  of  the  shroud  which  was  two  Inches.  The 
ambient  pressure  is  a  constant  and  is  taken  to  be  thirty  six 
pounds  per  square  foot.  The  change  in  thrust  is  found  to  be  a 
loss  of  0.5  pounds  force  from  a  nozzle  thrust  of  14.7  pounds 
force.  This  is  a  3.7  percent  loss.  If  the  ambient  pressure  is 
chosen  to  be  zero,  however,  a  thrust  gain  of  0.4  pounds  force  is 
calculated,  which  is  a  2.2  percent  gain.  Thus,  only  under 
certain  ambient  conditions  does  the  shroud  increase  the  thrust 
over  the  nozzle  alone. 

The  run  time  on  the  Cray  XMP  for  eighty  thousand  time  steps 
was  2373  seconds.  The  efficiency  parameter  is  then  1.4  x  10”* 

Couette  Flow 

The  Couette  flow  problem  has  been  used  previously  as  a  model 
problem  to  test  codes  (4,18).  A  Couette  flow  consists  of  a 
stationary  lower  plate  and  a  moving  upper  plate  with  a  fluid 
between  the  plates.  For  the  purposes  of  this  work  the  fluid  Is 
air  and  the  flow  is  considered  two-dimensional  and  steady  state. 
The  characteristic  length  is  the  distance  between  the  plates.  The 
exact  solution  to  this  problem  for  an  incompressible  flow  is,  for 
the  velocity  distribution  given  by  Eq  39  and  for  the  temperature 
distribution  is  given  by  Eq  40.  The  pressure  and  the  coefficient 
of  viscosity  were  assumed  constant  and  the  density  is  calculated 
from  the  equation  of  state. 

Figure  17  shows  the  results  from  the  PANS  and  FANS  codes  for 
the  Couette  flow  problem.  Both  codes  are  in  very  good  agreement 
with  the  exact  solution.  Besides  providing  a  dry  run  for  the 
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Plguce  17  Couette  flow:  velocity  distribution 
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codes,  this  exercise  also  demonstrated  one  very  important  point. 
In  both  codes  the  flow  fields  were  initialized  so  that  forward 
differences  could  be  used  on  the  "split"  terms  and,  as  mentioned 
earlier,  this  should  relax  the  step  size  restriction.  Eq  49 
applied  to  this  problem  yields 

(AX*Li,  >  i.6 

The  step  size  used,  however,  was 

Ax*  a  OJ 

In  short,  initializing  the  flow  field  and  using  forward 
differences  on  the  split  terms  significantly  reduces  the  size  of 
the  streamwise  step  size. 

Supersonic  Boundary  Layer 

As  a  second,  more  severe  test  for  the  PANS  and  FANS  codes 
the  supersonic,  laminar  boundary  layer  flow  previously  worked  by 
Lawrence  and  Tannehill  (10)  will  be  solved. 

This  is  a  flat  plate,  zero  pressure  gradient  flow.  The 
freestream  conditions  were  the  following: 

M,  =  2.0 

Tj  =  399.0  Rankine 
Re^  =  1.65  x  10® 


The  characteristic  length  was  taken  to  be  one  meter  measured  from 


the  leading  edge  o£  the  plate. 

The  flow  was  initialized  at  x*  =  0.305  using  data  from  the 

boundary  layer  code  of  Cebeci  (6).  The  initial  data  was  then 

* 

marched  to  X  =  0.915  and  compared  to  data  at  the  same  location 
from  Cebecl's  code.  Step  sizes  used  were  ax*  =  0.001  and  Ay*  = 
0.1524  x  10"*  with  six  hundred  and  twelve  grid  points  in  the  x- 
direction  and  forty  grid  points  in  the  y-direction. 

Figure  18  presents  the  tangential  velocity  profile  computed 
with  a  single  sweep  of  the  PANS  and  FANS  codes  compared  to  the 
boundary  layer  code.  Figure  19  presents  the  temperature  profile. 
Clearly,  the  PANS  code  provides  a  more  accurate  solution  than  the 
FANS  code  for  a  single  sweep.  This  is  expected,  since  for 
a  single  sweep  the  pressure  vector  was  set  to  zero  in  the  PANS 
code  and  the  negative  flux  vector  was  set  to  zero  in  the  FANS 
code.  Thus,  the  PANS  code  only  neglects  a  portion  of  the 
upstream  traveling  information,  but  the  FANS  code  completely 
eliminates  this  information. 

Next,  the  FANS  code  was  applied  to  the  same  problem,  but 
using  global  iterations.  Figures  20  and  21  compare  the  FANS  code 
run  for  four  sweeps  to  the  PANS  code  run  for  one  sweep  and  the 
boundary  layer  code.  Now  that  the  minus  flux  term  has  been 
included  in  the  FANS  code,  the  solution  is  just  as  accurate  as 
the  PANS  solution. 

An  indication  of  the  relative  computer  effort  required  by 
the  two  codes  is  given  by  a  comparison  of  CPU  times.  The  figure 
of  merit  here  is  the  CPU  time  in  seconds  divided  by  the  number  of 
sweeps  divided  by  the  number  of  grid  points.  For  the  FANS  code 
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Figure  21  Boundary  Layer:  Temperature  profile  (four  sweeps) 
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this  was  4.6  x  10"4  and  for  the  PANS  code  4.7  x  1CT4 


Shock /Boundary  Layer  Interaction  (FANS ) 

As  another  test  the  FANS  code  was  applied  to  the 
shock/boundary  layer  interaction  problem  previously  solved  using 
the  explicit  code. 

To  keep  the  number  of  grid  points  to  a  minumum  stretching 
in  the  normal  direction  was  used.  The  stretching  was  determined 
using  a  geometric  progression  (5:8) 


AyB’  =  M**iB  r“" 


(68) 


where  r  is  called  the  common  ratio.  It  was  found  that  the 
algorithm  could  not  handle  a  large  amount  of  stretching.  The 
maximum  value  of  the  common  ratio  was  found  to  be  approximately 
1.05.  By  specifying  the  location  of  the  outer  boundary  and  the 
number  of  points  in  the  normal  direction,  the  minimum  y  step  size 
could  be  calculated  from 


=  JB*(l-r)/(l  -  r-) 


Using  the  same  domain  size  used  by  Beam  and  Warming  gives 


min  ~  2.9097  X  10- 


With  a  grid  using  four  hundred  and  eighty  five  points  in  the 
streamwise  direction  and  one  hundred  and  one  points  in  the  normal 
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direction,  solution  accuracy  would  be  of  the  same  order  as  Beam 


I  and  Warming's  solution. 


Figures  22  and  23  present  the  results  of  the  computation. 

As  can  be  seen  in  figure  22,  the  pressure  contour  plot,  the  code 
I  did  an  excellent  job  of  capturing  the  leading  edge  shock,  the 

i 

•  impinging  shock  and  the  reflected  shock. 

For  comparison  the  results  of  the  explicit  code  and  Thomas' 
results  have  been  plotted  on  figure  23.  Here,  there  is  quite  a 
J  difference.  The  region  over  which  the  pressure  rises  from  before 
the  shock  interaction  to  after  the  shock  interaction  is  smaller 
in  the  FANS  results.  There  is  also  a  "hump"  in  the  FANS 
,  pressure  and  a  smaller  hump  in  Thomas'  data,  but  the  exploit 
pressure  curve  levels  out.  Both  of  these  differences  are  attri¬ 
butable  to  the  damping  in  each  code.  The  explicit  code  contains 
artificial  damping  terms  and  Thomas'  code  contains  damping  which 
■j  results  from  the  upwind  differences,  but  the  FANS  code  does  not 
contain  added  damping.  This  means  that  the  FANS  code  will  cap¬ 
ture  the  shocks  with  less  smearing,  but  overshoots,  as  seen  in 
the  pressure  hump,  may  occur. 

.* 

The  code  took  sixty  eight  seconds  and  three  global  sweeps  to 
solve  the  problem  on  the  Cray  XMP .  This  gives  an  efficiency 
parameter  of  4.6  X  10”4  .  Compared  to  the  results  for  the  same 

problem  for  the  explicit  code,  the  FANS  code  is  14  times  faster. 
The  efficiency  parameter,  however,  is  about  an  order  of  magnitude 
larger  which  is  to  be  expected  since  the  FANS  code  was  not 
vectorized. 
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Appendix  A:  Dimensions  of  the  Shroud 


The  shroud  used  in  this  work  is  a  model  of  nozzle  block  3B 
used  in  the  experiments  conducted  by  Moran  (11).  The  block 
consisted  of  a  single  converging-diverging,  two-dimensional 
nozzle  cut  out  of  aluminum.  The  sides  of  the  shroud  were  formed 
from  plexiglass  and  the  top  and  bottom  wall  of  the  shroud  was 
made  from  removable  wood  blocks.  The  figure  on  this  page 
shows  the  shroud  dimensions  and  the  drawing  on  the  next  page 
the  original  nozzle  design. 
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Appendix  B:  Jacoblans 


In  this  appendix  the  elements  of  the  Jacobians  defined  by 
Eqs  55,  56,  59  and  60  are  presented.  These  Jacobians  are  the 
partial  derivative  of  a  vector  with  respect  to  a  vector  and  are 
therefore  matrices.  Each  element  will  be  denoted  by  the  symbol 
given  to  the  corresponding  Jacobian.  Subscripts  on  the  symbols 
denote  the  row  and  column  of  the  element. 

For  the  [Q]  matrix: 

Qu  =  0 

Qii  =  l 


Qu  =  o 
Qu  =  o 

Qa.  =  (1/2X(7-3)u,9hK7-1)»',3 

Qtt  =  £2 — C7~l)Ju* 

Qt»  =  -(7-l>* 

0f4  *  (7-1) 

Q«  =  -ttV 

Qm»/ 

Qu-u* 

Qm  *  o 
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041  *  4(7Bi 7/0  - 

0<*  =  (7Bi*//>)  -  (7-lX3u**-*-»**)/2 
Ou  *  -{7-1)u*t* 

Q44  =  7U* 

For  the  [QJ*-  matrix  the  following  are  defined: 

c,  =  (dm'/dU,1)  =  [t<7-1)/2»‘J  HE.7/.^  * 

0,  =  (a.Vau.1  »  Mt-0/2»‘]  [-u*//»*] 
c,  =  (a»‘/flu47  »  Mr- iV2»“]  [-*7/>“l 

c4  =  (a»7atO  =  w7-t)/2a7  [i/p7 
b  =  ((-r-i)/-r]u*  *  2a‘/r 
B,*  =  /.•«-((l/2XM1-l)f 

D  =  {KT-lJu'-VfAT*-!)}  .  »n/2 


the  elements  are 

Oil*  =  (1/4)  [(u**/ s*Hp*u*V 

Qw"  »  (1/4)  [(2u7»*Hp*u*,/***)ci-2^p*c1] 


Q»*  =  (1/4)  [/>*-{p*u*a/a**)]c1 
Qm*  =  (1/4)  [A(p  YV2)] C4 


QaT 

*  BQu*  ■ 

-  (b,7p)  H7-i)u7p 

-2cJ 

Qaa* 

=  BQa" 

♦  (Kh)  [(7-1X1 //>•)  *  SoJ 

Qa»" 

=  BQM* 

*  (2B,77)c. 

Qu 

=  bqm* 

*  (2E,7t)°, 

o»r 

=  Q11V 

-  B,!///.-) 

Qsa" 

=  QiaV 

Qm" 

=  Qi»V 

-  BiW) 

Q34" 

=  QM  V 

QY  *  Q|.+D  Ei+{[-(7-1)W p*TH4a*-*-2(7-l)u*)c1-{2(7-l)a*u7 />*)]/(t*-1)  -  /*//>*} 


04t+  =  Qia*D  ®i*{[(7— l)2^*/ />*  4a*Cj  2(7-1)11*02  2(7-lJa*//»*]/(79-l)} 

=  Qu*D  -  E,+(l4a*Cj  -  2(7-1)u*c8]/(72-1)  -  y’/p*} 

Q44"  =  Qh+D  -  E1+{[4a,-2(7-l)u*]c4/(72-l)} 


For  the  elements  of  the  [R]  matrix 


rs  'r^  VT?T  T.T  V  V  FM  FTT^rT^r  V.  ■."  R*  R»  m"  V  w  >■ »  en  wyw?  w?  w*j  y.^r.-yjyjyjr.rT.' 


lyir.’V’vnr 'w  >r  v 


i 


i 


Ru  =  o 

Ru  =  0 
Ru  =  l 
Rl4  =  0 
R*  =  -uV 

R»  =  ▼* 

R*»  =  u* 

Rf4  -  o 

R*|  =  [(y-l)u*7-*(7-3)r**]/2 

R Si  =  “(7-l)u* 

R*S  s  (3— 7}t* 

R*4  =  7-1 

R41  *  [-(7RtVp*H7-iX«**^*a)>* 

Ru  =  -{7-1)11  V 

R4S  =  (7*V)  -  (7-1Xu^3y**)/2 
Ru  =  7y‘ 

For  the  elements  of  the  [W]  matrix 

Wu  s  0 

W|J  5  0 

W||  SC  0 
Wh  as  0 
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s 


3 

» 

r. 


Vs  V\v 


i 


w«  =  -/./Re  [d^/pl] 
w«  *  M/Re  [0# //>•)] 

WM  =  0 

?  wM  =  0 

w«  =  -<4/3X/./Re)  a^7/>1 
W„  =  0 

w»  =  (4/3)U/Re)  dWA 

Wj4  =  0 

W„  -  -/./Re  [a^u-//.*)  -  MWA/A  -  (T/Prjajflp’Ar-DpT  -  Iu***t**/2^*1>J 
w„  =  (/./ReXl-r/Pr)  a^u7/>7 
W„  =  (M/ReX4/3  -  7/Pr)  a^**//) 

W„  =  WReXT/Pr)  «,»//>*) 


When  the  derivative  approximations  of  Eqs  63  and  64  are 
applied  to  the  left-hand-side  of  Eq  62,  a  block  tri-diagonal 
matrix  results.  A  block  tri-diagonal  matrix  is  a  matrix  which 
has  only  three  bands  of  non-zero  block  elements.  These  bands 
consist  of  the  diagonal,  sub-diagonal  and  supra-diagonal  block 
elements.  Each  block  element  is  a  matrix  which,  in  this  case, 
are  four  row  by  four  column  matrices.  Now,  let  [A]  be  a  diagonal 
block  element,  [B]  a  sub-diagonal  block  element  and  [C]  a 
supra-diagonal  block  element.  The  elements  are  then  given  by 

[A] j  =  [Qlj1  -  [(Ax*/2Re(Ay*)*)  Mi 

[B] j  =  -{(ax72aj*XR]h  -  (Ax72ReUy*)*XM  *-jOMm] 

[C] j  =  (Ax*/2Ay*XR3j*i  -  (Ax72Ee(Ay*)a)U*rMiLiXw]j+i 
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Block  19 

'The  purpose  of  this  work  was  to  obtain  the  solution  to  three 
supersonic  flow  problems  using  three  different  numerical  techniques 

First,  a  shock/boundary  layer  problem  is  solved  using 
MacCormack's  explicit  technique.  Then,  using  the  same  technique 
a  shrouded  rocket  nozzle  problem  is  solved.  These  two  problems 
showed  that  the  explicit  scheme  required  many  minutes  of  computer 
time  to  solve. 


In  order  to  explore  more  efficient  codes  to  solve  the  shroud 
problem,  space  marching  algorithms  were  studied.  A  space 
marching  algorithm  using  flux-splitting  in  the  streamwise 
direction  was  applied  to  an  approximate  form  of  the  Navier-Stokes 
equation.  FI ux -sp 1 i t t i ng  combined  with  a  global  iteration 
approach  should  allow  the  shroud  problem  to  be  solved  with  a 
space  marching  algorithm.  The  flux-splitting  code  was  applied  to 
two  supersonic  flow  problems  and  very  good  results  were  obtained. 


